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abstract 


The  Transit  Improvement  Program  (TIP)  satellite  clock  is 
controlled  by  the  IPS  (Incrementally  Programmable  Synthesizer) 
and  flight  computer  subsystems.  Together  they  provide  a synthesis 
of  frequency  offset  and  frequency  drift  that  can  be  used  to  com- 
pensate for  such  errors  in  the  satellite  crystal  oscillator.  To 
do  this,  the  ground  software  must  estimate  the  oscillator  offset 
and  drift  and  then  compute  the  proper  control  parameters  to  be  in- 
jected into  the  satellite  to  steer  the  clock.  The  system  provides 
a resolution  of  control  of  1 part  in  10^3  in  frequency  and  1 part 
in  1013  per  day  in  drift.  Estimation  of  oscillator  offset  and 
drift  from  high-resolution  pseudo-random  noise  epoch  measurements 
is  accomplished  by  a discrete  Kalman  filter,  based  upon  a three- 
state  model  with  continuous  random  walks  in  frequency  and  fre- 
quency drift  as  the  driving  noise  in  the  oscillator.  A ground 
software  program  has  been  provided  to  implement  the  Kalman  filter 
and  compute  the  control  parameters  required  to  steer  the  satellite 
clock  to  the  reference  ground  clock.  The  programs  are  written  in 
Fortran  IV.  Complete  listings  of  the  software  and  operating  pro- 
cedures are  provided. 
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1.  INTRODUCTION 


The  Incrementally  Programmable  Synthesizer  (IPS)  subsystem 
on  the  Transit  Improvement  Program  (TIP)  satellite  series  provides 
a new  method  of  satellite  clock  control,  heretofore  untried  with 
quartz  crystal  oscillators.  The  subsystem  synthesizes  a frequency 
offset,  the  magnitude  of  which  is  progranmable  by  digital  input. 

The  hardware  is  described  in  Ref.  1. 

The  IPS  is  driven  by  the  high-quality  quartz  oscillator 
(5  MHz) , and  it  outputs  an  offset  frequency  that  is  controlled  by 
two  digital  registers  (A  and  B).  If  f()  is  the  oscillator  fre- 
quency, the  output  frequency  is 

f * fo  [‘-Will  • <» 

Thus  by  manipulation  of  the  A and  B registers,  the  output  fre- 
quency can  be  controlled  directly. 

The  hardware  clock  on  the  old  navigation  satellites  (OSCAR's) 
is  maintained  by  counting  the  oscillator  frequency.  The  traditional 
method  of  epoch  control  is  to  selectively  delete  or  add  counts  to 
compensate  for  oscillator  frequency  variations.  This  is  usually 
called  the  "clock  delete  system."  The  counter  on  the  TIP  satellites 
is  driven  by  the  IPS  output  frequency  rather  than  by  the  oscillator 
frequency.  With  this  system,  epoch  control  can  be  maintained  by 
direct  frequency  control  using  the  A and  B registers. 

By  making  small  adjustments  in  A,  the  IPS  output  frequency 
can  be  set  slightly  high  or  low  to  "steer"  the  epoch  error  to  zero 
between  settings.  Also,  by  allowing  the  A register  to  change 
continually,  the  crystal  aging  drift  and  flicker  noise  can  be  com- 
pensated for.  It  may  be  seen  from  Eq.  (1)  that 


f - A . (2) 

A B 


Ref.  1.  L.  Rueger  and  A.  G.  Bates,  "Frequency  Synthesizer 
for  Normalizing  the  Frequency  and  Time  Scale  of  Crystal  Clocks," 
28th  Annual  Symposium  on  Frequency  Control,  1974,  pp.  395-400. 
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This  relationship  can  be  used  effectively  to  compensate  for  oscil- 
lator drift  because  of  the  presence  on  TIP  of  a general-purpose 
flight  computer. 

The  IPS  hardware  is  interfaced  with  the  flight  computer  so 
that  the  computer  can  exercise  direct  control  on  the  A and  B 
registers  in  real  time.  With  this  system,  very-high-resolution 
clock  control  can  be  achieved.  The  IPS  hardware  as  implemented 
has  an  inherent  resolution  of  frequency  control  as  large  as  1 part 
in  10ll,  set  by  the  size  of  the  A register  (14  bits)  and  the  mag- 
nitude of  the  B register  contents.  By  using  the  flight  computer 
program  to  manage  the  IPS  registers  in  real  time,  the  resolution 
of  frequency  control  is  1 part  in  10^,  and  the  drift  control 
resolution  is  1 part  in  10-^  per  day.  The  flight  computer  program 
for  IPS  control  is  described  in  Appendix  A. 

Another  TIP  satellite  subsystem  provides  the  capability  for 
high- resolution  clock  epoch  transfers  to  the  ground.  This  is  the 
pseudorandom  noise  (PRN)  time  pulse  modulation.  The  PRN  code,  con- 
sisting of  a special  fast  phase-modulation  pattern  either  2*2  or 
2*5  chips  long,  is  transmitted  in  synchronism  with  the  satellite 
U.T.  clock.  By  using  dual-frequency  transmission  to  correct  for 
first-order  ionospheric  effects,  the  PRN  modulation  provides  the 
capability  for  epoch  recovery  with  a precision  of  several  nanosec- 
onds. The  accuracy  of  recovering  the  mean  satellite  clock  epoch 
for  a single  combined  doppler-PRN  pass  is  limited  to  about  20  ns 
by  the  uncertainty  of  the  satellite  position. 

The  complete  closed-loop  system  for  clock  control  using  the 
IPS  is  shown  in  Fig.  1.  In  this  typical  feedback  system  observa- 
tions are  made  of  the  epoch  error  and  controls  are  applied  to  null 
the  error  signal  to  zero.  The  interesting  features  of  the  system 
are  that  the  loop  is  closed  through  the  ground  software,  and  the 
time  constants  for  applying  controls  are  rather  slow  compared  to 
most  real-time  feedback  loops.  For  practical  reasons  that  have  to 
do  with  the  satellite  being  in  view  of  ground  stations,  the  con- 
trols can  be  applied  at  most  once  per  day  with  the  current  Transit 
ground  system.  Also,  for  the  same  reason,  epoch  measurements  can 
be  received  only  about  four  times  per  day.  Therefore  this  loop  is 
quite  slow  by  most  standards. 

It  is  assumed  that  the  epoch  measurements,  6tc,  input  to  the 
filter  and  control  program  are  not  tne  raw  PRN  measurements.  They 
must  be  processed  to  the  extent  of  providing  a single-parameter, 
mean-epoch  error  at  closest  approach.  The  parameter  is  estimated 
by  simultaneously  navigating  the  receiver  position  and  determining 
the  clock  error.  The  navigated  (rather  than  true)  receiver  posi- 
tion la  used  to  determine  propagation  time,  since  this  removes  the 
main  satellite  ephemeris  errors. 
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If  PRN  epoch  measurements  are  not  available,  any  clock  mea- 
surements will  suffice.  Also,  satellite  frequency  estimates  from 
tracking  runs  can  be  input  to  the  filter  program  to  supplement  the 
clock  measurements.  A Kalman  filter  is  used  to  estimate  the  state 
of  the  satellite  oscillator  from  the  epoch  and/or  frequency  mea- 
surements. It  is  convenient  to  set  up  the  calculation  of  the  IPS 
control  parameters  as  a recursive  estimation,  where  the  program 
is  run  with  the  latest  accumulated  data  just  prior  to  the  time 
scheduled  for  injecting  the  controls. 

A good  general  description  of  all  the  PDP-11  ground  software 
for  communicating  the  satellite  is  given  in  Ref.  2.  The  flight 
computer  software  system  is  described  in  Ref.  3.  The  remaining 
sections  of  this  report  describe  the  IPS  filter  and  control  pro- 
gram. 


Ref.  2.  C.  Marvin,  "An  Introduction  to  TIP-II  Satellite 
Digital  Operations,"  APL/JHU  SDO-4318,  January  1976. 

Ref.  3.  J.  M.  Whisnant  and  R.  E.  Jenkins,  "Post  Launch  and 
Operational  Flight  Computer  Programs,"  APL/JHU  SDO-4268,  Vols.  I 
and  II,  May  1976. 
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2.  THE  KALMAN  FILTER 


Many  of  the  results  in  this  section  are  standard  results 
from  linear  filtering  theory.  They  are  included  for  completeness 
and  to  aid  persons  not  familiar  with  the  Kalman  filter. 

In  the  filter,  we  model  the  IPS-driven  clock  output  as  be- 
ing determined  by  three  state  variables:  epoch  error  (fit),  IPS 

output  frequency,  and  frequency  drift.  The  satellite  hardware 
clock  counts  2*3  cycles  of  the  IPS  output  frequency  divided  by  12. 
When  the  IPS  output  frequency  is  exactly  nominal  (exactly 
4 999  577.60  Hz)  the  counter  underflows  precisely  6103  times  in 
two  minutes. 

We  define  the  last  two  states  to  be  relative  frequency  off- 
set from  nominal  and  relative  frequency  drift: 


fif 


f 


actual  nom 
f 

nom 


actual 

f 

nom 


f - 4 999  577.60  Hz 
nom 


(3) 


Then  the  propagation  equations  of  the  state  variables  as  a func- 
tion of  elapsed  time,  t,  are 


• (Ti ' T0r 

fit  - fit  + fif  (t.  - T_ ) + f 

T1  To  To  1 0 To  2 


fif  “fif  + f (T,  - T_ ) 

'l  To  To  1 0 


and 


i - f 
Ti  To 


(4) 

(5) 


(6) 
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The  transition  matrix  for  the  state  vector 


We  assume  that  the  above  system  Is  driven  by  Integrated 
white  Gaussian  noise  (i.e.,  random  walks)  in  both  frequency  and 
frequency  drift,  so  that  the  state  propagation  equation  for  small 
time  intervals  is  given  by 


(10) 


(ID 
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are  the  time  integrals  of  the  white  driving  noise  in  frequency  and 
frequency  drift,  respectively,  over  the  time  interval  - tq- 
Equation  (9)  is  an  approximate  imbedding  of  our  continuous  system 
in  a discrete  time  step  model.  This  is  valid  as  long  as  the  time 
step  is  infinitesimally  short.  Equation  (9)  has  been  included  to 
make  it  intuitively  clearer  how  the  noise  is  assumed  to  be  driving 
the  system.  When  - if)  becomes  larger,  the  \>2  and  V3  noises  mix 
among  the  states  in  a more  complex  manner. 

Since  our  time  steps  will  be  finite,  we  need  a better  ap- 
proximation than  that  afforded  by  Eq.  (9).  We  get  this  by  inte- 
grating Eq.  (9)  over  a span  of  infinitesimal  time  steps.  When  this 
is  done,  and  using  the  transition  property  of  <J>,  we  find  that 


Equation  (12)  is  usually  called  the  "matrix  superposition  integral," 
and  the  second  term  is  the  accumulated  noise  from  Tq  to  t^.  Equa- 
tion (12)  describes  a forced,  linear,  first-order  system  in  which 
the  dynamics  are  modeled  as: 


x (t)  - F(t)  x (t)  + I v2(t) 

lV3(t) 


(13a) 


where  F(r)  is  related  to  iKt.tq)  by 


— 4>(t,Tq)  - F (t)  4>(t,Tq) 


(13b) 


The  covariance  of  the  white  driving  noise  is 


q(*tu)  K { l v2(A)  J [0  v2(p)  v3(p)]>.  6(A  - p)  | 0 o 
w3(x) 


(14) 


- 13  - 


TNi  JOHNS  MCSWNS  UNIWlNSrtv 
APPLIED  PHYSICS  LABORATORY 

LAUNCL  MWVUM) 


where  6(X  - p)  Is  a Dirac  delta  function.  The  covariance  matrix 
of  the  accumulated  noise  after  a long  time  step  is 


QtVV 


dX 


T1 

J dp  (K^.X)  q(X,y)  ^(t^p) 


(15) 


An  explicit  evaluation  of  Eq.  (14)  with  At  ■ t^  - t^  gives: 


q(At)  - 


021dlil+(J2j[AT^  02IAt1 


2 3 


3 20  2 2 3 8 3 6 


»22  V 


, 2 lilL 

3 6 


2 (AX)  + 0 * iAlL-  0 2 lAti 


3 3 


2 IAt^ 
3 2 


3 2 


Oj  (At) 


(16) 


Equations  (12)  and  (16)  describe  the  modeled  IPS  output 
state  during  the  time  spans  between  chan  es  in  the  control  parame- 
ters A,  B,  and  K.*  When  these  parameters  change,  the  frequency  and 
frequency  drift  undergo  sudden  changes  not  described  by  the  state 
transition  matrix.  This  does  not  represent  a serious  problem. 

Given  a state  estimate  at  a time  just  before  a change  in  the  con- 
trol parameters,  we  can  obtain  an  estimate  of  the  state  immediately 
following  the  change,  as  we  will  see  in  Section  3.  In  tfye  process, 
of  course,  there  will  be  significant  changes  in  $f*  and  f*,  the 
estimated  frequency  and  frequency  drift.  There  will  also  be  changes 
in  the  covariance  of  the  state  estimate,  I.  However  the  changes  in 
the  covariance  matrix  can  be  shown  to  be  negligible  and  they  will  be 
Ignored  in  the  filter. 


In  running  the  filter,  we  assume  that  we  have  some  estimate 
of  the  state  at  a time  tj,  following  the  last  change  in  control 
parameters,  and  based  upon  all  previous  measurements: 


+K  is  defined  in  Appendix  A. 
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(17) 


Further,  the  covariance  matrix  of  the  error  In  this  estimate  is 
also  generated  at  time 


We  are  given  a number  of  measurements  of  the  epoch  error 
and/or  frequency  offset  at  times  t^,  and  we  must  produce  an  esti- 
mate of  the  output  state  at  some  future  time,  t, , when  new  con- 
trol parameters  are  to  be  injected  into  the  satellite. 

The  measurements  are  corrupted  by  noise  that  we  assume  is 
uncorrelated  from  measurement  to  measurement,  so  that  the  measured 
quantities  are 


z “ H x + v(x  ) (19) 

L i 


where  for  a time  error  measurement, 


H - (1  0 0) 


(20) 


and  for  a frequency  measurement, 


H - (0  1 0) 


(21) 


The  covariance  "matrix"  of  the  scalar  noise,  v(t^),  is  the  variance 


Rfv^)]  - c6t2 


(22) 
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for  a time  measurement,  or 


R[v(Tt)]  - Of2 


(23) 


for  a frequency  offset  measurement. 

The  optimal  solution  of  this  estimation  problem  Is  provided 
by  a Kalman  filter,  as  described  In  Ref.  4.  The  filter  starts  with 
estimates  of  the  state  and  of  the  covariance  matrix  at  some  time 
tl,  given  all  the  previous  measurements  up  to  some  time  tp;  that 
Is,  we  know 


X*  (tl|tp) 


and 


E<tJV 


Given  a measurement  at  time  t^,  we  first  propagate  these 
estimates  to  the  time  of  the  measurement: 


x*  (tJV  " ♦(t1  ~ tl)  x*  <tJtp)  (24) 

E(ti|tp)  - *(t1  - tl)  £(tJtp)  ♦+(t1  - tl)  + Q(r±  - tl)  . (25) 

We  then  compute  the  measurement  residual 

ZCrjTp)  - Z - H x*  (rjTp)  (26) 

and  the  optimal  filter  gain 

. * — 1 

K(rt)  - r(T1|tp)  H [R  + H T(t1|tp)  H+]  (27) 

Ref.  4.  R.  E.  Kalman,  "New  Methods  In  Wiener  Filtering," 
Proceeding*  of  First  Symposium  on  Engineering  Applications  of  Ran- 
dom Function  Theory  and  Proabllltv.  J.  Wiley,  1963. 
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and  finally  update  the  state  estimate  and  covariance  matrix 

X*  Ct±  | T±)  - x*  (t±  | Tp)  + K(Tt)  Z(t1|tp)  (28) 

Z(T1IT1)  - ^(TjTp)  - K(Tl)  H Z(t1|tp)  . (29) 


This  process  Is  continued  for  each  measurement.  We  note 
that  the  matrix  Inversion  in  Eq.  (27)  is  trivial  in  our  case  since 
R and  H £ IT*"  are  scalars.  Note  also  that  the  values  for  H and 
R used  to  process  a given  measurement  depend  on  the  measurement 
type  (fit  or  fif). 

After  the  measurements  have  been  processed,  the  state  esti- 
mate is  propagated  up  to  the  next  injection  time,  t.,  using  equa- 
tions similar  to  (24)  and  (25).  To  do  this,  the  next  injection 
time  must  be  known;  hence,  the  data  should  not  be  processed  until 
the  IPS  injection  has  been  scheduled  so  that  tj  can  be  input  to 
the  program. 
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3.  COMPUTATION  OF  CONTROL  PARAMETERS 


Having  passed  the  data  through  the  Kalman  filter,  we  have 
an  estimate  of  the  clock  state  at  the  next  scheduled  Injection 
time,  Tj.  To  compute  the  required  control  parameters  for  that  In- 
jection, we  require  estimates  of  the  frequency  and  frequency  drift 
of  the  oscillator.  These  are  obtained  from  the  clock  state  esti- 
mate at  Tj  along  with  the  values  of  A,  B,  and  K at  Tj. 

We  start  from  the  values  of  these  parameters  at  tl,  the 
time  of  our  previous  Injection.  (Notice  that  tj  from  the  last 
program  run  becomes  tl  for  the  current  run.)  The  values  are  given 
In  their  "external  form,"  that  Is,  the  actual  numbers  that  are 
transmitted  to  the  flight  computer.  In  effect,  these  parameters 
are  then  converted  Into  an  "Internal  form,"  and  It  Is  this  Inter- 
nal form  that  occurs  In  the  control  equations  (1)  and  (2).^  The 
transformation  from  external  to  Internal  form  Is  given  by 


A -*■  A + 1 

B -*■  B + 1 

K -*>  K • 2 


(30) 

(31) 

(32) 


Since  K and  B do  not  change  with  time  once  they  are  in- 
jected, their  Internal  values  at  time  r,  are  maintained  until  Tj. 
However,  A does  change  with  time  as  shown  In  Eq.  (2),  and  its 
value  must  be  propagated  up  to  time  Tj.  The  actual  time  variation 
in  A Is  controlled  by  the  flight  computer  (see  Appendix  A,  Eqs. 
(A-l)  to  (A-3)).  In  the  flight  computer,  A is  incremented  every 
At  by  a value 


AA  - K 


(33) 


where  At  “ 128*120/6103  *2.52  a.  Hence,  the  difference  equation 
satisfied  by  A with  2.52-s  intervals  is 


tThla  is  an  Idiosyncrasy  of  the  hardware.  To  avoid  confusion  we 
have  treated  the  Internal  form  inside  the  software  so  that  it  can 
be  Ignored  by  persons  running  the  program. 
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n+1 


A + 
n 


— A 2 

2 n 


(128) 


(34) 


We  have  been  unable  to  solve  this  equation  exactly.  How- 
ever, a similar  difference  equation, 


A " A + 
n+1  n 


(128) 


2 An  An+1 


(35) 


should  be  an  adequate  approximation  to  (34)  for  these  purposes. 
This  is  exactly  solvable  as 


A 

n 


1 - 


(128) 


2 V 


(36) 


so  that  our  propagation  equation  for  A is 


A(Tj) 


A<V 


K«6103 
(120) (128) 


3 A(TL)(TI 


V 


(37) 


The  time  rate  of  change  of  A Is  manifest  as  a drift  In  the 
output  frequency  so  that  (see  Appendix  A) 


i f0  A f0  6103 

f - B J m — 2 3 

B A BA  (120) (128) 


(38) 


Therefore,  having  the  values  of  the  control  parameters  at  hand,  we 
can  arrive  at  an  estimate  of  the  oscillator  state  by  using  Eqs.  (1) 
and  (38): 


(39) 
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0 


1 - T (1  + T)]  (f 


f*  K (6103)  \ 
B (120) (128) 3 / 


(40) 


Values  of  A,  B,  and  K must  now  be  chosen  to  correct  the  out- 
put to  zero  frequency  drift  and  to  provide  an  output  frequency,  f8, 
that  will  steer  the  epoch  error  back  to  zero  over  some  time,  T: 


f - f 
s nom 

f 

nom 


(41) 


To  do  this,  we  require  values  of  A and  B so  that 


f 

s 


fo*  t1  - i (1  ♦ J>] 


(42) 


where  A and  B are  constrained  to  lie  in  the  ranges  (see  Ref.  1) 


6301  « A £ 13  470 


(43) 


6669  < B < 12  428 


(44) 


The  algorithm  for  choosing  these  numbers  is  as  follows.  First, 
compute  the  control  constant: 


C 


1 

B 


(l+i) 


(45) 


Next,  find  a (noninteger)  value  of  B so  that  C obtains  when  A 
has  the  value  6301: 


B “ (l  - cX1  + 6301  ) * <46> 
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I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 


Choose  as  candidate  values  of  B the  integers  straddling  this 
value,  computing  A from 


(47) 


where  [B]  is  the  integer  either  greater  than  or  less  than  B. 
Finally,  choose  the  A,  B pair  that  fits  within  the  constraints  of 
(43)  and  (44). 

The  value  of  K is  then  chosen  to  compensate  for  the  esti- 
mated drift: 


(48) 


Finally,  the  output  state  estimate  is  updated  to  the  cor- 
rected values  at  time 


(49) 


--61-0-3-  ,«0 

(120) (128) 


(50) 


The  approximate  equalities  in  Eqs.  (49)  and  (50)  are  not 
exact  since  the  control  parameters  have  finite  precision. 

The  external  form  of  the  new  control  parameter  can  be  found 


by  the  transformation 

A > A - 1 

(51) 

B B - 1 

(52) 

31 

K -*■  K • 2ja  . 

(53) 
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4.  RUNNING  THE  SOFTWARE  PROGRAM 


The  software  implementation  of  the  algorithms,  written  in 
Fortran  IV,  consists  of  a main  program  and  nine  subroutines.  The 
program  listings  are  shown  in  Appendix  B. 

The  subroutine  FILT  implements  the  Kalman  filter.  The  sub- 
routine PRED  propagates  the  state  estimate  up  to  the  input  injec- 
tion time.  These  routines  call  matrix  multiplication  routines 
MM3  and  MMT3.  Calculation  of  the  control  parameters  is  done  by 
subroutine  IPS.  Card  input  and  output  are  performed  by  the  rou- 
tines INPUT,  TASK,  and  OUTPUT,  which  were  isolated  to  simplify 
changes.  For  the  same  reason,  a separate  BLOCK  DATA  routine  is 
used  to  set  the  noise  parameters  assumed  by  the  filter. 

The  program  is  designed  to  run  on  the  IBM  360-40,  and  in- 
puts are  via  punched  cards  in  the  present  version.  One  group  of 
input  cards  (the  "state"  cards)  is  punched  out  by  the  program  each 
time  it  is  run,  and  these  cards  become  the  input  state  cards  for 
the  next  run.  The  remaining  input  cards,  containing  the  measure- 
ment data,  must  be  newly  prepared  before  each  run. 

The  use  of  cards  is  one  area  that  the  Naval  Astronautics 
Group  (NAVASTROGRU)  may  wish  to  change.  It  may  be  more  convenient 
to  use  one  or  more  disk  files  as  the  input  and  output  stream,  where 
the  state  output  card  images  are  written  by  the  program  onto  the 
disk  to  be  read  by  the  subsequent  run.  Also,  the  data  on  some  of 
the  measurement  cards  could  come  directly  from  the  program  that 
processes  the  satellite  clock  data.  It  may  eventually  prove  con- 
venient to  have  this  program  write  a disk  file  for  the  measurement 
card  images  that  would  be  read  by  the  IPS  control  program.  In 
order  to  effect  this  type  of  change,  all  that  is  required  is  to 
change  the  subroutines  INPUT,  TASK,  or  OUTPUT  to  read  or  write  a 
different  tape  number,  and  define  the  appropriate  files  in  the  JCL 
setup. 


The  first  cards  in  the  input  deck  are  the  six  state  cards 
that  describe  the  previous  state  estimate  at  time  t^.  The  formats 
of  the  six  cards  are  as  follows: 

CARD  1 FORMAT  (115,  F15.8)  - gives  the  last  injection  time, 

tl 

IDAYL  - day  number 

SECL  - U.T.  seconds  of  day 
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CARD  2 FORMAT  (3D20.13)  - gives  clock  state  at 

DELT  - epoch  error  <microseconds) 

F - IPS  frequency  (ppm  offset  from  5 MHz) 

FD  - frequency  drift  (ppm  per  day) 

CARDS  3,  FORMAT  (3D20.13)  - gives  covariance  matrix  at 

4,  P - 3 * 3 covariance  matrix  (internal  units) 

5 

CARD  6 FORMAT  (D20.13,  2110)  - gives  control  parameters  at 


A,  B,  K - control  parameters  (external  form) 

With  the  exception  of  the  covariance  matrix,  all  variables 
on  the  state  cards  are  in  "external"  units.  Internally,  epoch 
error  is  converted  to  seconds,  frequency  to  fractional  offset  from 
the  nominal  frequency,  frequency  drift  to  fractional  offset  per 
second,  and  A,  B,  and  K to  their  internal  forms  as  per  Eqs.  (30), 
(31),  and  (32).  The  covariance  matrix  is  input  in  these  internal 
units. 


Following  the  state  cards  are  the  series  of  task  (measure- 
ment) cards  with  the  following  format: 

TASK  CARD  FORMAT  (2115,  2F15.8)  - one  for  each  measurement 
MEAS  - flag  indicating  type  of  measurement 
IDAYM  - day  of  measurement 
SECM  - measurement  of  time  of  day  (seconds) 

Z - value  of  measurement 

If  MEAS  ■ 1,  the  card  represents  a time  error  measurement  at  epoch 
IDAYM,  SECM;  and  Z is  in  microseconds.  If  MEAS  ■ 2,  the  card 
represents  a frequency  measurement,  and  Z is  in  ppm  offset  from 
5 MHz.  Following  all  the  measurements  (which  must  be  time  ordered) 
is  a final  card  in  which  MEAS  ■ 0,  indicating  that  a control  in- 
jection is  to  take  place  at  time  IDAYM,  SECM.  (The  Z value  is 
ignored  on  this  card  and  may  be  blank. ) 

The  card  output  from  the  program  is  the  set  of  six  state 
cards  to  be  used  as  input  the  next  time  it  is  run.  The  printed 
output  is  generated  by  the  main  routine  and  the  IPS  subroutine. 
Figure  2 is  an  example  of  the  output. 
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IP5  control  Porr.HAM 


I NI T I * L VALUFS: 

l A ST  COMPUTE)  PFOrw: 

DAY*  27  ^FCONnS=  *3200. 000000 

eSTIMATF0  STATE  VAPIABLcS: 

time  cprdPs  -0.3*9  MICROSECONDS 

FREQUENCY*  -A*. *799901 73  PARTS  PER  MILLION  OFFSET 

f REOHE  NC Y DR  I CT  * -0.000000032  OAPTS  PCR  MILLION  OFFSET  PER  DAY 

rrjVARlANCF  MATRIX  l.NTERNAL  UNITS): 

0.999*33567792*0-13  0.276512621*1080-17  0.5879507*0*2050-23 
0.276512621*1080-17  0.92976585719320-2?  0. 200778206333 8D-27 
0.5879507*0*2050-23  0.20077820633380-27  0.31131*07355160-32 
IPS  CONTROL  PARAMETERS: 

A*  6505. 73911 18930000  8=  11838  K = 569 


T IMF  MEASUREMENT  AT  Ep0CH: 

0 AY  * 28  SECONDS*  0.0 

TIME  error*  -0.*23  MICROSECONDS 

FILTER  GAINS: 

K l = 0.998*323*156270  00  K2=  0.15333692910970-0*  K3  = 0.305**6022*7720-10 

FILTERED  state  E S T I MA TF : 

TIMF  ERROP*  -0.*23  MICROSECONDS 

FREQUENT  Y*  -8*. *79990159  PARTS  PER  MILLION  OFFSFT 

FREOUcNCY  ORIFT=  -0.000000027  PARTS  PER  MILLION  Offset  pep.  DAY 


INJECTION  AT  FPOCH: 

DAY*  28  HOl>o*  12  MINUTES*  0 SFCONDS*  0.0 

STATp  FSTIMATF: 

TIME  ERROR*  0.002  MICROSECONDS 

EREOUENCY*  -8*. *79990173  PARTS  PER  MILLION  OFFSFT 

FRFQUFNCY  DRIFT*  -0.000000027  PARTS  PER  MILLION  OFFSET  PER  0 AY 

VALUES  OF  IPS  PARAMETERS  AT  INJECTION  TIMF: 

A*  6529.32905086*0810  8*  11838  X*  569 

OSCILLATOR  STATE  AT  INJECTION  TIMF: 

FREQUENCY*  -0.000*622*1  PARTS  PEP  MILLION  OFFSET 

FREQUENCY  DRIFT.  -0.0000*6925  PARTS  PER  MILLION  OFFSET  PER  OAY 

STEERING  FOR  0.0  MICROSFCONOS  TIME  FRROR 

AT  86*00.000  SECONDS  PAST  INJECTION  TIME 

DESIRED  OUTPUT  FREQUENCY: 

FREQIJFNCY*  -8*. *80000022  PARTS  PER  MILLION  OFFSFT 

POSSIBLE  VALUF*  O'  IPS  PARAMETERS: 

A*  652*. 360222**63300  8*  11838  K*  569 

A*  *205. *919357869050  B*  11839  K*  569 


Fig.  2 Sample  Output  of  IPS  Program 
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CHOICE  OF  IPS  PAR  AMCTCRS  S 

A*  6 524. 3602  2 2446  3300  8 = 11838  K = 569 


STATF  ESTIMATE  AFTER  INJFCTION: 

LAST  COMPUTEO  FPOCH: 

OAY-  28  SECONDS*  43200.000000 

ESTIMATED  STATE  VARIABLES: 

TINE  ERROR-  0.002  MICROSECONDS 

FREOIIENC  Y = -84.480000022  PARTS  PER  MILLION  OF F SF T 

FREOUENCY  DRIFT*  -0.000000027  PARTS  PER  MILLION  OFFSET  PER  DAY 

COVARIANCE  MATRIX  (INTERNAL  UNITS): 

0.99943525897720-13  0.27651318979830-17  0.58795903161010-23 
0.27651318979830-17  0.9297677674777O-22  0.20078099290620-27 
0.58795903161010-23  0.20078099290620-27  0.31131814338360-32 
IPS  CONTROL  PARAMETERS: 

A*  6524.3602224463300  8*  11838  K*  569 


4-4 

Fig.  2 (cont'd) 
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The  BLOCK  DATA  subroutine  contains  seven  parameters  that 
may  need  to  be  changed  occasionally.  They  are: 


Program 

Name 

Analytic 

Notation 

(Eq.) 

Definition 

S2MD 

of2  (23) 

Variance  of  frequency  measurement  error 
(fractional  offset,  squared) 

S2M 

°st2  (22> 

Variance  of  time  measurement  error  (sec- 
onds, squared) 

Q2 

o22  (14) 

Variance  of  frequency  driving  noise  (frac- 
tional frequency  offset,  squared,  over  a 
1-s  interval) 

Q3 

o32  (14) 

Variance  of  frequency  drift  driving  noise 
(fractional  offset  per  second,  squared, 
over  a 1-s  interval) 

FNOM 

f (3) 

non 

Nominal  output  frequency  (ppm  offset  from 
5 MHz) 

STEER 

A nonzero  time  error  to  which  the  output 
will  be  steered,  if  desired;  zero  other- 
wise (seconds) 

TAU 

T (41) 

Time  in  which  output  time  error  is  to  be 
steered  to  STEER  (seconds) 

The  compiled  values  for  the  listed  parameters  are  based  on 
experimental  runs,  some  with  simulated  data  and  some  with  a real 
oscillator.  Some  of  the  values  will  depend  on  how  the  program 
is  being  run,  and  the  BLOCK  DATA  routine  will  have  to  be  recom- 
piled accordingly.  The  recommended  values  for  the  parameters  are 
described  below. 

With  PRN  recovered  epoch  errors,  the  value  of  S2M  is 
9 * 10~-6,  Using  this  value,  the  satellite  would  be  injected 
every  day,  and  the  program  would  be  run  once  per  day  with  a value 
of  TAU  of  86  400. 

With  epoch  errors  recovered  from  the  satellite  fiducial 
time  in  the  navigation  message,  the  value  of  S2M  is  4 * 10“^. 

In  this  mode,  the  program  will  need  to  be  run  less  often,  possibly 
once  per  week.  The  required  frequency  of  injection  will  be  left 
open  to  be  operationally  determined  by  NAVASTROGRU.  In  any  case. 
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the  value  of  TAU  should  be  compiled  accordingly.  The  actual  mea- 
surements entered  into  the  program  should  be  the  processed*  clock 
data,  averaged  over  each  pass.  The  corresponding  measurement  time 
should  then  be  the  time  of  closest  approach  for  the  pass. 

Bear  in  mind  that  the  program  runs  recursively  and  chrono- 
logically in  time.  This  means  that  the  program  should  always  be 
run  with  the  new  Injection  time  later  than  the  last  measurement 
time  in  the  input  deck.  Also,  if  there  are  no  measurements,  a new 
injection  time  must  still  be  supplied. 

If  a run  has  to  be  repeated  for  some  reason  (such  as  the 
scheduled  injection  not  being  made),  the  run  must  be  restarted  with 
the  old  state  cards,  not  with  the  state  cards  from  the  invalid  run. 
This  is  extremely  important  because  the  program  assumes  that  the 
IPS  parameters  on  the  state  cards  have  actually  been  injected  into 
the  satellite  at  This  is  the  only  way  it  can  predict  the  cor- 

rect clock  behavior  to  compare  to  the  measurements.  Therefore  to 
recover  from  any  problem,  the  program  run  sequence  must  be  re- 
started with  the  state  cards  for  the  last  valid  injection. 

If  no  state  cards  are  available  because  they  have  been  lost, 
or  if  the  satellite  clock  is  reset  completely  because  of  a flight 
computer  failure,  the  state  cards  will  have  to  be  generated  again 
to  restart  the  process.  This  is  done  by  selecting  any  convenient 
epoch  and  supplying  the  best  estimate  of  the  clock  state  and  its 
covariance  matrix  at  that  epoch.  For  the  frequency  offset  and 
drift,  these  values  can  be  obtained  by  extrapolating  from  the  last 
known  history  as  long  as  the  oscillator  has  not  lost  power.  If  the 
oscillator  has  lost  power  and  is  restarted,  then  a satellite  orbit 
fit  will  be  needed  to  obtain  the  estimates. 

In  the  absence  of  better  information  to  start  from  the  be- 
ginning, the  covariance  matrix  can  be  input  as  a diagonal  matrix: 


where 

-3 

* expected  error  in  clock  epoch  estimate  (10  if  clock 
was  just  reset) 


tCorrected  for  time  delays,  etc. 
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t~  * expected  error  in  offset  (10  ^ if  there  is  an  offset 
from  orbit  fit) 

e.  ■ expected  error  in  drift  (10  ^ if  there  is  drift  from 
3 orbit  fit). 

Once  the  program  is  run  with  valid  data,  the  covariance  matrix  will 
eventually  stabilize  to  a range  of  values  determined  by  the  measure- 
ment errors  and  the  assumed  filter  noise  parameters. 
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Appendix  A 

FLIGHT  COMPUTER  IPS  CONTROL+ 


By  the  nature  of  ita  design,  the  TIP  IPS  subsystem  cannot 
achieve  the  frequency  or  frequency-drift  resolution  necessary  for 
high-precision  time  control  with  fixed  settings  in  the  IPS  regis- 
ters. This  resolution  problem  has  nothing  to  do  with  the  stability, 
phase  drift,  or  temperature  sensitivity  of  the  circuitry.  Bench 
tests  have  shown  that  the  hardware  performs  excellently.  The  reso- 
lution is  limited  only  by  the  lengths  of  the  registers. 

Better  resolution  is  achieved  by  using  the  flight  computer 
to  manipulate  the  IPS  registers  in  real  time.  Without  a great 
deal  of  trouble,  a relative  frequency  resolution  of  10”^^  and  a 
relative  drift  resolution  of  10“13  per  day  can  be  obtained.  This 
is  sufficient  to  maintain  10-ns  accuracy  with  one  setting  a day. 
Whether  that  accuracy  can  be  realized  in  orbit  depends  largely  on 
the  performance  of  the  oscillator  and  the  precision  of  the  PRN 
clock  epoch  error  measurements. 

Preliminary  tests  with  the  engineering  model  of  the  TIP  os- 
cillator showed  a real  potential  for  better  than  100  ns  control. 

To  exploit  the  full  potential  of  the  system,  a flight  computer  pro- 
gram to  control  the  IPS  registers  has  been  written  to  produce  a 
resolution  of  10“^^,  The  program  will  be  loaded  with  new  IPS  con- 
trol parameters  periodically  (once  a day)  via  the  ground  station 
computer.  The  parameters  will  be  generated  from  epoch  or  frequency 
measurements  on  the  ground  using  the  IPS  filter  and  control  program. 

Four  values  are  input:  A,  B,  AF,  and  K.  A and  B are 
16-bit  integers.  A goes  directly  to  the  IPS  A register;  B goes 
directly  to  the  B register.  AF  and  K are  32-bit  integers. 

AF  is  the  fractional  part  of  A and  represents  the  fraction 
of  time  that  the  A register  should  hold  (A  + 1) . 

Effective  IPS  A - [3  • @ 


tThis  appendix  la  based  on  an  internal  memorandum  by  R.  E.  Jenkins, 
"Flight  Computer  IPS  Control,"  APL/JHU  S1A-88-74,  16  December  1974. 
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Out  of  every  128  tocks,  the  IPS  A register  should  hold  (A  + 1) 
for  .AF  x 128  tocks  and  should  hold  the  value  A for  128(1  - .AF) 
tocks.  We  pick  128  for  convenience  since  multiplication  by  128 
shifts  the  decimal  point  seven  places  in  AF.  This  means  we  pick 
up  the  first  seven  bits  from  AF  as  the  counter  for  changing  the  A 
register. 

K determines  the  secular  change  in  the  A register  to  com- 
pensate for  crystal  aging  drift.  Every  128  tocks,  the  value  of  AF 
is  permanently  incremented  (or  decremented)  by  AD: 

AD  - KA2  . 


Each  time  AD  is  added  to  AF  in  the  active  program,  the  carrys  from 
the  most  significant  bit  of  AF  increment  or  decrement  A by  one. 

The  sign  of  the  drift  correction  will  be  carried  by  K.  To 
simplify  the  program,  it  is  preferred  that  the  value  of  AD  does  not 
cause  a carry  in  a single  128-tock  interval  with  a reasonably  high 
drift  rate.  That  is,  we  want  the  secular  change  to  A to  be  less 
than  one  count  in  128  tocks,  and  let  the  accumulation  of  the 
changes  in  AF  carry  into  the  IPS  A register.  This  sets  the  maxi- 
mum correctable  drift  to  about  10~®  per  day,  which  is  more  than 
enough . 


Also,  we  want  the  precision  in  drift  correction  to  be  lO-^-^ 
per  day.  To  achieve  the  required  precision  in  calculating  AD  over 
the  full  range  of  values  of  A (6300  to  13  470)  and  B (6668  to 
12  427),  we  need  to  calculate  KA^  and  hold  AF  in  31-bit  precision. 
A practical  scaling  to  do  this  is 


"■*(158)  • <*-!> 


With  this  scaling,  K is  defined  as  follows. 

From  the  formula  for  the  IPS  transfer  function,  to  correct 
for  a relative  drift,  f/fQ,  the  required  increment  in  A per  At  is 


+We  have  defined  the  unit  of  time  in  the  flight  computer  clock  as  a 
"tock,"  equal  to  19.662461  ms.  The  flight  computer  maintains  a 
U.T.  clock  by  counting  tocks  generated  by  a computer  interrupt  com- 
ing from  a hardware  frequency  counter. 
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AD  - - f-  BA2  At 
0 


(A-2) 


For  At  - 128  tocks.  At  - 128  x 120/6103  - 2.516795019  s.  Thus 


f _ ,,  120 

AD  - - -j-  B<128) 


\128  / 


(A- 3) 


which  defines  K,  a number  considerably  less  than  one. 

To  transmit  K to  the  flight  computer,  which  has  16-bit 
words,  we  change  it  to  a double-bit  precision,  31-bit  binary  frac- 
tion: 


K - 


32  • 16  1 


Therefore  K-(f/f0)(B)  (128)3  120/6103  x 231  and  will  be  trans- 
mitted as  a right- adjusted  31-bit  integer.  The  range  of  this 
integer  will  be  from  1 to  about  210  000  for  the  worst-case  drift. 

A negative  value  for  K will  be  carried  as  a double-precision 
integer,  with  bit  32  indicating  sign. 

With  the  above  scaling,  the  value  of  (A/128)^  can  be  com- 
puted and  stored  in  single  precision  (16  bits),  since  the  maximum 
value  of  A is  16  272.  This  considerably  simplifies  the  calcula- 
tion of  AD  in  the  program. 

A detailed  description  of  the  IPS  management  program  and 
how  it  fits  into  the  overall  flight  computer  software  system  is 
given  in  Ref.  2.  Figure  A-l  is  a simplified  flow  chart  of  the  pro- 
gram. The  program  is  re- initialized  every  time  new  IPS  parameters 
are  injected  from  the  ground. 
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Frequency  counter  interrupt 


Fig.  A-1  Simplified  Flow  Chert  of  IPS  Program 
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IPS  CONTROL  PROGRAM 

PROGRAMMER:  4.  GOLDFINGER  JHU/APL 

DIMENSION  XI 31 .P( 3,3  ) 

DOUBLE  PRECISION  DEL T ,F ,FC. FC ,FNOM, X, P, Z, S2M, 02,03. STEER. TAU. A , 2 
DOUBLE  PRECISION  TL , TC . SE CL . SECM, DELD AY, TM,T I ,S2M0 
DOUBLE  PRECISION  SEC 

COMMON  70/  S2M,Q2,Q3,FNOM, STEER, TAU.S2MD 
10  FOR  4A  T I 20HI IPS  CONTROL  PROGRAM  I 
20  FORMAT ( I6H01NITI AL  VALUES:  I 
30  F0RMATI21H0L4ST  C0MPUIE3  EPOCH:  I 
AO  FORMA  T I 6M  DAY*  ,I3,UH  SECONOS*  .F12.6) 

Al  FGRMAT16H  DAY*  ,I3,BH  HOUR*  ,13.11H  MINUTES*  .13, 

1 l 1H  SECONDS*  ,F 12.61 
50  FORMAT127HOEST1MATED  STATE  VARIABLES:! 

60  F0RMAII13H  TIME  ERROR*  ,F10.3,13M  M ICROS ECONOS » 

70  FORMA  T ( 1 2H  FREQUENCY*  ,F1A.9,25H  PARTS  PER  MILLION  OFFSET! 

80  FORMAT  1 1 8H  FREQUENCY  DRIFT*  ,FtA.9, 

l 33H  PARTS  PEP  MILLION  OFFSET  PER  OAY ) 

01  FORMAT  < 36H  COVARIANCE  MATRIX  (INTERNAL  UNITS):  ) 

82  FORMATIIH  ,3020.13/1H  ,3020.13/lH  ,3020.13) 

83  FURMATI2AH  IPS  CONTROL  PARAMETERS:  I 

HA  FORMAT I AH  A*  ,F20.13,5H  8=  ,!6,5H  K*  ,16! 

90  F0RMATI20 HO  INJECT  ION  AT  EPOCH:) 

100  FOR"AT(27HOT IME  MEASUREMENT  AT  EPOCHS) 

110  FORMAT  1 32  HOF REQUE NCY  MEASUREMENT  AT  EPOCH!) 

120  FORMATI2SHOF ILTEREO  STATE  ESTIMATE:) 

130  FORMAT! 16H0STATE  ESTIMATE:) 

1 AO  FORMAT (32HOSTATE  ESTIMATE  AFTER  INJECTION:) 

150  FORMAT! 3 l HO* ******* ************ *•»******♦) 

READ  IN  DATA 

CALL  INPUTIICAYL,  SECL , OEL  T,  F , F I) ,P  , A,  1 6,K) 

WRITE (6,10) 

WRITEI6.150) 

WR 1TE (6,20) 

WRITE (6, 30) 

WRITE (6, AO)  IOAYL ,SECL 
WRITE (6,50) 

WR 1TE 1 6, 60)  DELT 
WRITE (6, 70)  F 
WRITE (6,80)  FO 
WRITE (6,81) 

WRITE (6,82)  P 
WRITE (6,83) 

WRITEI6.8AI  A, 19 , K 

CONVERT  DATA  TO  INTERNAL  FORM 

FNOM*  5»D0*( 1 .06 A FROM  ) 

X( l) *0ELT*1 .0-6 

X( 2) *5.00*11 .06*F l/FNOM-l.DO 

X( 3)*5.D0*F0/(86A00.Q0*FNCM) 

TL-SECL 
FC*X ( 2 ) 

TC*SLTL 
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FILTER  MEASUREMENTS 

200  CALL  TASK<MEAS,10AVM,SECN.2> 

MRITE  <6.1501 

IF (MGAS.EQ.O)  MR  I TE  t 6.901 
IF<MEAS.E0.1I  WR I TE ( 6. 100  I 
IF (ME  AS.ro. 2)  WRI TE ( 6. 1 10  I 
)F(MEAS.EO.OI  GO  TO  800 
WRITE  16.40)  IOAYM.SECM 
I F ( ME  AS.  E C.  1 I WR  I TE  ( 6 , 60  I 2 
IF (Me  AS .FQ. 2 I WR!TE(6.70>  2 
IFIMEAS.rQ.il  2=2*1. 0-6 
IF (ME  AS. EO. 2 I 2*5.00*(l.06*2) /FNOM-l.DO 
DELDAVIDAYM-1DAYL 
TM»ScCM*8640J.00*DEL0AY 
CALL  FILT(TL,TM,X,P,2,MEAS| 

WR ITE 16,1201 
DELT*X(1)*1.D6 

F * FNOM* ( 1 .00* A ( 2 1 1/5.00-1.06 
F 0*8 6 400.  DO*  FNl)M*X<  3 1/5.00 
WRITE (6,60)  OELT 
WRITE(6.70)  F 
WRITE  (6. <101  FO 
TL*TM 
GO  TO  200 

PREDICT  STATE  AT  TI 

800  f HR*SECM/3600.00 

SEC*S  ECM- 3600.00* IHR 
M I N*S  F.C/60.00 
SEC*SCC-60.00*MIN 
KR1TE(6,41)  IDAVM.IHR.MIN.SEC 
DE  LOA  Y* 1 0 AYM— I DAY  L 
T l *SE  CM*  86400.D0*DEL0AV 
CALL  PREDITL ,T I.X  »P  ) 

WRITE (6,1301 
OEL T»  X ( 11*1.06 
F»FN0K*(1.00*X(2I 1/5.00-1.06 
FO *86400. 00* FNOM* X< 3 1/5.00 
WR|TE(6,F0)  OELT 
WR ITE (6,70)  F 
WRITE (6, HO)  FD 

COMPUTE  IPS  PARAMETERS 

OELT*X ( 1 1 
F*XI  2 ) 

FO*X( 31 

CALL  IPS(TC.TI,OELT,F,FO.A,IB.K) 

XII ) *DELT 
X ( 2 ) * F 
XI3I-F0 
C 

C CONVERT  DATA  TO  EXTERNAL  FORM 

C 

DELT-Xt  11*1.06 

F*FN'JM*(  1 .D0»X(2I  >/5. 00-1.06 
F 0*86400. DO* FNOM* XI 3 1/5.00 
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OUTPUT  DATA 

WRITE (6, ISO) 

WRITE  (6. 160) 

WRITE (6,301 

WRITE (6. *0)  IDAYM.SECM 
WRITE (6.50) 

WRITE  16.60)  DEL  T 
WPITE (6,70)  F 
WRITE (6,80)  FO 
WR  ITt (6,81) 

WRITE(6,8?)  P 
WR  ITE  ( 6,  83) 

WRITEI6.8A)  A.I8.K 
WRITE (6,150) 

CALL  OUTPUT ( I OAYM ,SECM, CELT , F ,F0, P, A, IB,KI 

STOP 

END 
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BLOCK  DATA 

COMMON  / Cl/  S2M,02.03,FN0M,STEER.TAU.S2MD 
DOUBLE  PRECISION  S2M ,Q2 ,0 3, FNGM, STEER , T AU ,S2N0 
DATA  S2MD/ 1.0D-22/ 

DATA  S2M/ 4.0D-10/ 

TATA  02/ 2.00D-29/ 

DATA  03/  4.0D-40/ 

DATA  F NON/ -84 . 430 0/ 

DATA  STEER/O.OD*)/ 

DATA  TAU/ 864 00.00/ 

END 

SUBROUTINE  INPUT! 1DAYL»S53L»4ECT»6  60.7  1.9B  2' 
DIMENSION  P ( 3*31 

DOUBLE  PRECISION  DEL  T , F ,F  D,  P , A , SE CL 

INPUT  FCRMA1S  FOLLOW 

10  FCRMATII15,F15.8I 
2C  FORMA  T( 3D 20. 13) 

30  FORMAT  Iu20. 13*2110) 

READ  IN  OATA 

RE  AO (S. 101  ICAYL.SECL 
RE  AO ( 5 ,20 ) DEL  T,F  »FD 
DC  100  1*1.3 

100  READ!  5.20)  P ( I , 1 ) .P  II  ,2)  , PI  I ,31 
RE  AD  I 5 ,30 ) A.1H.K 
RE  TURN 
END 

SUBROUTINE  TASKIMEAS.IOAVM.SECM.Z) 

DCUHL E PRECISION  Z.SECM 

INPUT  FORMAT 

10  F0RMATI21 15.2F1S.B) 

READ  OATA 

RE  AD (5, 10)  ME  AS , I DAVM.SECM.Z 

RETURN 
ENO 

SUBROUTINE  OUTPUT ( I DAYP, SECM.DELT .F ,F0,P ,A , I B.K) 
DIMENSION  PI  3, 3) 

DOUBLE  PRECISION  DEL T ,F , F D, P, A , SECM 

OUTPUT  FORMATS  FOLLOW 

10  FORm.ATII  15, F 15.8) 

20  FORMAT(3D20. 13) 

30  FURMATID20. 13,21 101 

C 

C WRITE  OUT  OATA 

C 

WR IT  t I 7, 1 01  IDAYM.SECM 
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WRITE  IT ,20)  DELT,F,PD 
DO  100  I*lt3 

100  WRITE  I 7 » 20 ) P(!,l >,P(I,2I,P< 1,3) 
WRITE (7*30)  A.IB.K 
RETURN 
END 

SUBROUTINE  MM3<4,B,AB) 

DIMENSION  A(3,3),8(3,3),4B{3,3) 
OOUBLE  PRECISION  A,B,AB 
DO  100  1=1,3 
00  100  J=l,3 
AB ( I . JI=O.ODO 
DO  100  K=I,3 

100  A8I I , J!*AB( I . J>»AII ,K)*B(K, JI 

RETURN 

END 

SUhRUUTINE  MMT3( A,S( ABTt 

DIMENSION  Al 3,31 ,B! 3,3), ABTI 3,31 

OCUHLE  PRECISION  A, B , ABT 

00  100  1=1,3 

00  100  J= 1 , 3 

AHT ( l , JI =0.000 

DO  100  K=  1, 3 

100  ABTt I ,JI=ABTI I,J|«A( I,K)*B(J,K| 
RE  TURN 
END 
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SUBROUTINE  F ILT I T L,  T M,  X ,P  , Z , ME  AS  1 
COMMON  /C1/S2M.02  ,03  , FNCM ,ST E ER , T AU, S 2M0 

01  MEN  SI  ON  X(3).P(3,3I,E<3,3),PUI3,3>,PHII3,3),PX(3,3I,M3I 
DOUBLE  PRECISION  X, P , F , PU ,PH I ,P X, DENOM, S2M , 02 , 03 . K , ZT , Z ,DT,F , FO 
OCUOLE  PRECISION  FNO M , S TE ER . T AU , F C , TL , T M , S2M0, TOEL 
10  FORMAT ( ' OF  I L TER  GAINS:*) 

20  FORMAT!*  Kl  = *,020.13,*  K2  = *,020.13,*  K3*  *,020.131 

C 

C GET  TRANSITION  MATRIX 
C 

PH  1 1 l ,11  = 1.00 

PHI  ( t ,2>  = TM-TL 

PHI ( 1 • 3 ) =0. BOO* ITM-TL1**2 

PHI(2,1)=0. 000 

PHI (2 ,21=1.000 

PH II2,3)  = TM— TL 

PH  I ( 3 , 1 ) =0.0D0 

PHII3, 21=0. 000 

PH  1(3,31=1.000 

c 

C PROPAGATE  OLO  COVARIANCE  MATRIX 

C 

CALL  MMT3(P,PH1,PXI 
CALL  MM3 ( MHl ,PX,P) 

TDEL=  TM-TL 

P(l,l I=P ( l,ll*Q2*TQEL**3/ 3 .00*0 3* TOEL **5/20. DO 
P( l , 2 I * P ( 1 ,2 1 *02*  TO EL* *2/ 2. 00* 03* TO EL **4/ 8. DO 
PI  1,3  1 = P(  l, 3 l *03* TOEL** 3/ 6.00 
P I 2, 1 1 *P ( 2, 1 M02* TOEL**2/2.DO*03*TOCL**A/e.OO 
PI  2,2 ) = P(2,2I*O2*T0EL*O3*T0EL**3/3.OO 
P (2,3  I *P( 2,3 1*03* TOEL **2/ 2. 00 
PI3,l 1 *P( 3, 1 l*03*TDEL**3/6.00 
P( 3,2  1 =P I 3, 2 I *03* TOEL **2/ 2. 00 
PI3«3 )=P( 3,3 l*03*T0EL 
C 

C PROPAGATE  STATE 

C 

DT=X(  l)*PHI I 1,21* XI 21* PHI  (1 ,3)*X(31 
F = X 1 2 l*PH  I ( 2 , 3 I *X 1 3 ) 

FD«X( 3) 

XI II=0T 

X I 2 1 * F 

X I 3 1 = FD 
C 

C COMPUTE  F ILTER  GA  INS 

C 

IFIMEAS.EO.il  GO  TO  250 
IF  IMF AS.EQ.2 ) GO  TO  2/0 
250  0EN0M=S2M*P( 1,11 
K I l)=P(l,l)/0EN1M 
RI2)»P12,11/0EN0M 
K(  3)  = PI  3,  ll/OENJM 
GO  TO  300 

270  0ENf'M  = S2M0*PI2,21 
K ( 1 1 * P C 1 , 2 1 / OENOM 
K(21*P(2,2 1 /OENOM 
K|  3 I = P I 3, 2 1 /DENOM 
300  CONTINUE 
C 
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GET  MEASUREMENT  RESIDUAL 

IF (ME AS«EQa 1 1 GO  TO  330 
I F (ME AS.EQ.  2 1 GO  TO  370 
3 30  ZT-Z-XCl) 

GO  TO  400 
370  ZT-Z-XC2) 

GET  NEW  STATE  ESTIMATE 

40C  X(l)*Xlll*K(l)*ZT 
XI2)-XC2)»M2)*ZT 
XC3)-X(3I*«(3I*ZT 
WRITE  1 6 » 10) 

WRITE(6,20I  K 

UPDATE  COVARIANCE  MATRIX 

DO  600  1-1,3 
DO  600  J» 1,3 
600  EC  I, J 1-0.000 
EC  1,1  1-1.000 
EC2, 21-1.000 
EC3. 31-1. 000 
IF1MEAS.E0.1)  GO  TO  630 
IFCMEAS.E0.2I  GO  TO  670 
630  Ell,l  1-El  1,1  l-KC 1 I 
E C 2,  l ) ®E l 2, 1 I-KC2  I 
E I 3,  l )-E ( 3, 1 )-K( 3 I 
GO  TO  700 

670  EC 1,21-FC 1,2 l-KC  1 ) 

EC2.2 )-E ( 2«2  l-K C 2 I 
EC  3,2 l-EC  3,2 )-M 3 I 
700  CALL  MM3CE,P,PU) 

00  800  1-1,3 
DO  800  J*  l,  3 
800  PC  I.JI-PUCI, Jl 
RETURN 
ENO 
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SUBROUTINE  PREDi TM, T I , X,P > 

COMMON  /Cl/  S2M, 02,03, FNOM,STEER.TAU,S2MD 
DIMENSION  X{ 3), XI 13) , PI  3,3), PX( 3, 3), PHI (3. 3) 
OOUBIE  PRECISION  X, X I ,P ,P X, PH  I , TM ,T I , S2MD, TDEL 
DOUBLE  PRECISION  S2M ,02 ,03, FNCM, STEER. TAU. PC 
C 

C GET  TRANSITION  MATRIX 

C 

PH  1(1.11*1.00 
PHI(1,2)rTI-TM 
PHI  11 • 3) «0. 5D0* IT l-T  M)P*2 
PH  1 1 2 • 1 )*  0.00 
PHI (2 .2 1 ■ 1.00 
PHI(2,3)»TI-TM 
PHI (3,1 I-O.DO 
PHI ( 3 , 2) *0.00 
PH  1(3.31*  1 .00 
C 

C PROPAGATE  STATE 

C 

XI(l)«X(l)+PHl(l,2)*X(2)»PHI(l,3)*X(3> 

XI (2) *X(2)»PHI { 2 , 3) *X( 3 ) 

XI ( 3 ) *X( 3 ) 

X ( 1 I * X I ( 1 > 

X ( 2) * X I ( 2 ) 

X ( 3 ) * XI I 3 ) 

C 

C PROPAGATE  TRANSITION  MATRIX 

C 

CALL  MMT3 (P.PHI.PX) 

CALL  MM3( PHI ,PX,P) 

TDEL*  T I-TM 

PCl.l )*P(1,1 ) ♦Q2*TOEL**3/ 3,00*Q3*TDEL**5/20.00 

P( 1, 2 ) »P ( l»2)*Q2*TOE  L**2/2. D0*Q3*TDEL**A/8.D0 

PI1.3 l*P( 1.3l»03*T0rL**3/6.D0 

P( 2. I ) *P( 2, 1 )♦ U2 * TO EL** 2/ 2.00*03* TDEL  *•4/8.00 

P(2,2)*P(2,2)*02*T0£L*03*T0EL**3/3.D0 

P(2,3l»P(2,3)»03*TDEL**2/2.00 

P( 3,1 )*P( 3, l I *03* T0EL**3/6.D0 

P ( 3.2 l«P( 3.2 )»03*TDEL**2/2.00 

P(3,3)*P(3,3)»U3*TOEL 

RE  TURN 

END 
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SUBPOUT  INC  IPS! TO. T I ,DELT ,F,F0, A, IB.KI 
DOUBLE  PRECISION  OELT ,F ,FD,FC, FNOM, S2M,Q2,Q3, STEER, TAU 
OOU9L E PRECISION  F0«  FOG, FOP «F00P*CIPS ,B, fcK, EK AL T . A. AALT , ENINT . POUT 
OCUBL  T PRECISION  T0,TI,TINT,EKINT,FCIUTP,A0,S2M0 
COMMON  /Cl/  S2M,Q2,03,FN0M, STEER, TAU, S2MD 
10  FORMAT  I A4.H0VALUES  OF  IPS  PARAMETERS  AT  INJECTION  TINE* I 
15  FORMAT (36H0QSCI LL  ATCR  STATE  AT  INJECTION  TIME*) 

20  FO RMA T ( 12H  FREQUENCY-  .F1A.9.25H  PARTS  PER  MILLION  OFFSET) 

30  FORMAT ( 18H  FRFOUhNCY  DRIFT-  ,Fl*.9, 

l 33M  PARTS  PER  MILLION  OFFSET  PER  DAY  I 

31  FORMAT! 1 A HOST  EES  I NG  FOR  ,F10.3,2*H  MICROSECONDS  TIME  ERROR  ) 

32  FORMA T I AH  AT  ,F11.3,28M  SECONDS  PAST  INJECTION  TIME) 

35  FORMA  T ( 26H0DE  S I Rr  0 OUTPUT  FREQUENCY  1 1 

50  FORMA  T(35HOPOSSIBLE  VALUES  OF  IPS  PARAMETERS!) 

60  FOHMATIAH  A-  .F20.13.SH  A-  .16. 5H  K-  ,16) 

70  FORMAT  (26H-JN0  A*H  SOLUTIONS  IN  RANGE  I 
8C  eORMA  T I26H0CM0ICE  OF  IPS  PARAMETERS*! 

C 

C CONVERT  DATA  TO  INTERNAL  FORM 
C 

A-Atl.ODO 
IB-IBFl 
0-19 
EK-K 

EK-EK *2.000**1-31  I 
PROPAGATE  A TO  T1 
TINT-TI-TO 

ENINT -TINT/2. 51679501800 
A-A/l  1 .ODO-EN INT*EK* A/ 163 BA. 0001 
WRITE (6,101 
A0-A-1.000 
I BO- 1 0-1 
KO-K 

WRITE  16,601  AO, I BO, KO 

GET  OSCILLATOR  FREQUENCY  AT  TI 

Cl  PS- 1 .000-1 1.000*1 .000/ A l/B 
FO-(FfI.OO»/CIPS-1.00 
FOP-F  NOM* 1 1 . 00»F0 )/ 5.00-1 .06 
WRITE  16, 151 
WRITE  16,20)  FOP 

GET  OSCILLATOR  ORIPT 

FDO-I FD-I F0*1 .00) *EK  *6103 .00/ 1 B* 120. DO* 126.00*6 3 ) )/C I PS 
FOOP- 116*00.  DO*FNOM*FOO/5.  00 
WRITE <6, 30)  FOOP 

GET  OESIREO  OUTPUT  FREQUENCY 
NOTE!  OESIREO  FOO-O 

STEERP-STFEP*I.0E6 
WR  I Y t (6,311  STtEFP 
WRITE  16,32)  TAU 
TOUT -ISlfFR-OfLTI/TAU 
FCUTP-FNUM*! l.OJ* FOUT 1/5. 000-1.06 


44 


F 


I 


r 


i 


i 

i 

i 


i; 


TMt  JOHNS  HOPKINS  UNIVtBSITV 

APPLIED  PHYSICS  LABORATORY 

LAUREL  MARYLAND 


WR  ITE  (6,351 
WRITE(6,20>  FOUTP 
C 

C GET  IPS  PARAMETERS 
C 

C INITIAL  A AND  B SO  THAT  F-FOUTS 

c 

Cl  PS*  IFOUTM  .001/  (F0*1.D0) 

B*l.D0/(l  .OO-CIPSI 
iB*fuo/63oi.no 
A« 6/  I IB-BI 
I BAL  T * I H*  1 
AALT>B/( IBALT-BI 
C 

C GET  K sn  THAT  F00UT»0 

C 

E K *-F CO* I B* l 28.00 **3*120. CO/ 1 I FO* 1.00 1*61 03. 001 *C I PS 
EKALT=-FDU*IBALT* 128.00**3*120.00/1 (F G»1 .00 1*6103 .001 *C IPS 

c 

C CHOOSE  A,B,X 

c 

MR ITE (6.501 
AOA-l.DO 

i eo= i b- i 
KO*EK*2.DO**31 
WR|TE(6,60>  AO.IBC.KO 
Afl*AA  LT-  1 .00 
I BO* I BALT-1 
KO*FK ALT*2.D0**3l 
WR ITE (6.601  AC.IBQ.RO 
1*1 
I X-l 

IF  UA.LT.6301.).0R.( A.GT.134  70. M 1-0 
IFUAALT.LT.  6301.  !.0R.(  AALT.CT.13A70.  It  IX*0 
I F ( ( I 6. LT. 6665 ». OR. ( I B. GT . 12 428 1 1 I >0 
IFUIFALT.lt  .6669  I.Ofc  . ( I b AL T .GT  . 1 2628  1)  IX-0 
1FUI  .EO.OI  .ANO.I  IX. EC. 01  I GC  TO  200 
IF(( l .EO.OI.ANO.I IX.NE.OII  GO  TO  300 
IFUI  .NF.OI  .ANO.I  IX.EO.O)  ) GO  TO  400 
IFU.LT.AALT  I ICH-1 
IFIA.GE.AALT ) !CH*2 
GO  TO  500 
200  WR ITE (6,701 
GG  TO  600 
300  ICH*2 

GO  TO  500 
400  I CH» 1 

GO  TO  500 

500  IF(ICH.EQ.l)  GO  TO  550 
A*AALT 
16*18  ALT 
EK*EX  ALT 
550  WR ITE 16, 80) 

AO* A- 1 . 00 
I BC*1 B-l 
KIJ»EX*2.0D0**31 
WR ITE (6,601  AO.IBO.KO 

c 

C GET  STATE  AFTER  INJECTION 

c 
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•■IB 

C IPS*  1 .00-1 1 .00* 1 .DO/AI /• 

EK-K0*2.000**(-3l I 
F«l*0»l.D0l*C IPS- 1.00 

FO-FDO*C  IPSMFUM  .00  t•EK•«10^.00/(••1^0.00•i2•.00••^l 

C CONVERT  DATA  TO  EXTERNAL  FORM 
C 

60C  AkA-1.00 
IB-IR-l 

EK«ER  *2.00**31 

K«EK 

RETURN 

ENO 


i 
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